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ABSTRACT 

At  re-entrant  corners  the  pressure  has  a  singularity  for  incompressible  viscous  flow.  In 
fluid  flow  computations  there  are  geometries  that  have  re-entrant  corners,  and  for  which 
it  is  needed  to  provide  an  appropriate  value  for  the  pressure  at  such  a  corner  when  a  finite 
difference  method  dealing  with  the  primitive  formulation  is  used.  In  this  paper  we  address 
the  problem  of  finding  an  efficient  strategy  for  computing  pressure  values  at  a  re-entrant 
corner  which  applied  to  Strikwerda’s  second-order  numerical  method  for  solving  the  Stokes 
and  Navier-Stokes  equations.  The  pressure  at  the  corner  is  regarded  as  a  double  valued 
function.  Also  we  examine  Moffatt’s  solution  for  the  Stokes'  problem  near  a  step  where  the 
pressure  becomes  unbounded  as  the  re-entrant  corner  is  approached.  We  show  that  this 
strategy  models  very  well  the  pressure  singularity  making  the  computation  more  amenable 
and  efficient. 
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NUMERICAL  TREATMENT  OF  THE  PRESSURE  SINGULARITY 
AT  A  RE-ENTRANT  CORNER 
Gerardo  A.  Ache’ 

1.  Introduction 

\ 

For  plane  steady  incompressible  viscous  flow,  it  is  known  that  in  geometries  with 
corners  the  flow  close  to  the  corner  will  be  dominated  by  the  Stokes’  solution  and  the 
governing  equation  for  the  stream  function  will  be  the  biharmonic  equation  [8]  i.e. 

V4rp  =  0  ,  (1.1a) 

with  V'  satisfying  the  non-slip  boundary  conditions 

dib 

^  =  ^°  ’  dn  =  °  '  (Ll6) 

As  pointed  out  in  [2],  [7]  and  [8],  corner  flows  lead  to  a  singularity  in  the  vorticity  and 
the  pressure.  As  examples  of  such  problems  having  geometries  with  corner,  that  one  sees 
often  in  applications,  are  flow  in  a  cavity  and  flow  in  a  channel  with  a  step  (or  a  pipe 
with  a  sudden  enlargement  of  the  cross  section).  For  flow  in  a  cavity,  in  most  of  the 
finite  difference  schemes,  the  singularity  for  the  pressure  or  the  vorticity  may  be  avoided 
if  it  is  assumed  that  the  presence  of  such  singularities  affects  the  solution  only  in  a  small 
neighborhood  of  the  corner.  We  find  a  different  situation  for  flow  near  a  step,  that  is  a 
re-entrant  corner.  Even  if  the  singularity,  in  the  pressure  or  vorticity,  is  assumed  to  be 
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negligible  it  is  needed  to  provide,  appropriate  values  for  the  pressure  or  vorticity  to  evaluate 
pressure  or  vorticity  derivatives  at  grid  points  near  the  corners.  An  excellent  survey  about 
vorticity  strategies  at  a  re-entrant  corner  can  be  found  in  [5]. 

In  this  paper,  we  address  the  problem  of  providing  appropriate  pressure  values,  at 
a  re-entrant  corner,  when  a  finite  difference  scheme  for  the  primitive  formulation  of  the 
Stokes  or  Navier-Stokes  equation  is  used.  As  an  example  of  such  a  scheme  we  examine 
[10].  We  show  that  using  the  approach  adopted  in  this  paper,  we  are  able  to  compute 
numerical  solutions  for  the  Stokes  equations,  in  domains  with  re-entrant  corners,  for  which 
the  numerical  pressure  accurately  models  the  pressure  singularity. 

2.  The  Fluid  Flow  Problem  Near  a  Sharp  Corner 

In  this  section  we  examine  Moffatt’s  solution  [8]  for  plane  flow  near  a  corner.  The 
stream  function  ip  is  defined  in  polar  coordinates  by 


where  if  the  corner  is  at  (x0,  y0)  then  z  =  z0  +  r  cos  8  and  y  =  y0  +  r  sin  0.  We  assume 
that  the  corner  is  given  by  —  a  <  0  <  a,  and  the  governing  equation,  in  polar  coordinates, 
is  given  by 


.  j  _  d*ip  2  dzip  1  d2ip  1  dip  2  d4ip 

*  ~  dr*  +  r  dr3"  ~  ^Ifr2  +  dr  +  H  dr2dO 2 

_2_dhp_  j.3V_0  (22) 

r3  drdB2  +  r*  d82  +  r4  d84  '  1  ’ 

In  [7]  Lugt  and  Schwiderski  have  shown  that  the  equation  (2.2)  admits  separable  solutions 


of  the  form, 


M’J)  =  r'  f\{8)  , 


(2.3) 


where 


fx(0)  =  A\ sin  X0  +  B\  cos  A 9  +  C\  sin(A  -  '2)6  +  D\  cos(A  -  2)0  ,  (2.4a) 

for  A  ^  1  and  , 

/  i(0)  =  A  i  sin  20  +  B\  cos  20  +  C\9  +  D\  .  (2.46) 

Thus  the  general  solution  to  (2.2)  can  be  written  in  a  series  of  the  form  , 

OO 

<6  =  <6o+£tfnrA»/A„(*)  ,  (2.5) 

n=  1 

with  the  eigenvalues  A„  satisfying 

0  ^  /2e(Ai)  ^  iEe(A2)  ^  ••• 

In  order  that  ip  satisfies  the  boundary  conditions  (1.16)  it  is  required  that  An  satis¬ 
fies  an  eigenvalue  equation.  For  symmetric  slow  motions  (B\  —  D\  =  0)  this  equation 
becomes,  (for  An  /  1) 

sin2Q:(An  —  1)  =  (An  -  1)  sin  2a  ,  (2.6a) 

while  for  skew  symmetric  motions  (.A*  =  C*  =  0)  , 

sin2a(An  -  1)  = -(A„  -  1)  sin2a  .  (2.66) 

The  pressure  P\  can  be  recovered  using  the  stream  function  ip  and  the  Stokes  equa¬ 
tions.  Neglecting  an  arbitrary  constant  one  arrives  at  (for  A  ^  1) 

P  \{r,9)  =  -4ArA-2C^  cos(A  -  2)0  ,  (2.7a) 


for  symmetric  motions,  and 


for  skew  symmetric  motions.  In  general  the  final  pressure  is  given  by  the  superposition  of 
symmetric  and  skew  symmetric  motions,  for  all  A  satisfying  (2.6a)  or  (2.6b).  For  flow  near 
a  step,  where  2a  =  37r/2,  we  can  combine  (2.6a),  (2.6b)  so  that, 


(An-l)(-l)"+'=Sm(A„-l)^ 


(2.8) 


where  for  n  odd  the  An’s  are  associated  to  symmetric  motions  and  for  n  even  with  skew 
symmetric  motions.  In  Table  I  we  list  the  first  four  eigenvalues  for  this  problem. 


TABLE  I 


Symmetric  and  Skew-symmetric  Eigenvalues 


n 

Eigenvalues 

wm 

1.54448 

CB 

1.90853 

KB 

2.62926+0.23125i 

Bfl 

3.30133-f0.31584i 

The  final  pressure  for  the  Stokesian  flow,  near  a  re-entrant  corner,  can  be  written  as 

p(r ,«)  =  ,  (2.9) 

n 

where  the  coefficients  C\n,D\n  in  P\n  are  determined  using  the  boundary  conditions.  One 
sees  from  (2.7a),  (2.9)  and  Table  I  that  the  pressure  becomes  unbounded  as  r  approaches 
the  re-entrant  corner  and  the  leading  term  is  approximately  of  the  order  0(r-0-456). 

For  the  fully  non-linear  two  dimensional  Navier-Stokes  equations  it  is  possible  to 
analyze  the  solution  near  a  corner  using  the  technique  described  in  [2).  This  approach 
is  based  on  a  semi-analytic  perturbation  analysis  which  is  valid  near  a  corner  for  small 


Reynolds  numbers.  The  use  of  this  technique  may  be  useful  to  study  the  nature  of  the 
singularity  of  the  pressure  near  a  re-entrant  corner  for  -small  values  of  R  . 

S.  Numerical  Treatment  of  the  Pressure  Singularity  at  a  Re-Entrant  Corner 

Although  the  pressure  at  a  re-entrant  corner  may  be  unbounded,  a  finite  value  is 
needed  for  the  finite  difference  scheme  at  such  a  corner  to  evaluate  the  pressure  gradient 
at  adjacent  grid  points.  A  similar  situation  is  encountered  when  the  stream  function- 
vorticity  formulation  is  used  since  the  vorticity  is  unbounded  at  the  re-entrant  corner  and 
yet  a  finite  value  is  needed  to  evaluate  first  and  second  derivatives  at  adjacent  grid  points. 
Some  corner  strategies  for  the  vorticity  when  a  finite  difference  approach  is  used  have  been 
discussed  in  [4],  (5],  [6],  [9]  . 

In  this  paper  we  provide  a  finite  difference  strategy  for  the  pressure  applied  to  the 
Stokes  and  Navier-Stokes  equations.  The  numerical  method  that  we  used  for  that  purpose 
is  based  on  a  second-order  finite  difference  scheme  to  solve  the  Stokes  and  Navier-Stokes 
equations  [10],  [11].  This  method  deals  with  primitive  variables  formulation,  i.e.,  it  uses  the 
velocity  and  the  pressure  as  dependent  variables.  Then  the  laplacian  and  convection  terms, 
given  in  conservative  form,  are  discretized  using  standard  centered  difference  formulas, 
while  the  pressure  gradient  and  the  continuity  equation  are  discretized  using  regularized 
centered  difference  formulas  [10],  e.g. 


r\  - 

~  «  ^yO«2  -  -h26y+62_u2  ,  (3. Id) 

where  (6io>£yo)»  (^x-,^-),  (6z+,<5y+)  are  centered,  backward,  and  forward  differences  in 
the  x  and  y  direction,  respectively.  Formulas  (3.1)  are  used  in  the  pressure  gradient  and 
the  continuity  equation.  To  determine  the  pressure  on  the  boundaries  cubic  interpolation 
formulas  are  used,  e.g.  along  a  boundary  given  by  t  or  j  equal  to  zero, 

Poj  =  3(pi,y  -  P2,j)  -  P3,y  ,  (3.2a) 


P.,o  =  3(pi,i  -  pi, 2)  ~  Pi, 3 


(3.26) 


To  solve  the  resulting  finite  difference  equations  an  extended  S.O.R.  iterative  proce¬ 
dure  [11]  was  used.  This  numerical  method  is  very  efficient  and  accurate.  It  has  been  used 
to  compute  flow  in  a  spinning  and  coning  cylinder  [12] . 

Since  the  numerical  scheme  requires  values  of  the  pressure  at  the  solid  wall  and  since 
these  values  are  computed  using  extrapolation  formulas,  we  may  compute  the  pressure  at 
the  corner  in  two  different  ways  using  such  extrapolation  formulas.  For  each  of  the  two 
directions  we  let  Pc,x  and  Pc,y  be  the  pressure  value  determined  by  extrapolation  in  the  x 
direction  and  y  direction.  Assuming  that  the  re-entrant  corner  is  located  at  (xo,yo),  we 
have 

P  c,z  :=  p(x 0,yo) 

=  3(p(x0  4-  hx,y0)  -  p(x0  +  2hI,y0))  +  p(x0  +  3hx,i/0)  ,  (3.3a) 

and 

pe,v  :=  p(*o,yo) 


=  3(p(x0,  y0  +  hy)  -  p{x 0,  y0  +  2 hy))  +  p(x0,  y0  1-  3/iy)  .  (3.36) 
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Here  hx,hy  are  the  grid  sizes  in  the  x  and  y  directions  respectively.  Notice  that  since  these 
values  of  Pc<x  and  Pc>y  do  not  need  to  be  equal  the  pressure  at  a  re-entrant  corner  may 
be  regarded  as  a  double  valued  function.  Therefore  to  evaluate  the  pressure  gradient  at 
adjacent  grid  points  we  use  Pzx  to  compute  the  approximation  to  dp/dx  at  (xo  +  hx,yo) 
and  Pc>y  to  compute  the  approximation  to  dp/dy  at  (xo,yo  +  /ij/).  This  strategy  works  well 
in  computation  at  small  and  moderate  Reynolds  numbers.  This  approach  applies  also  to 
axi-symmetric  and  three-dimensional  computations  which  include  flow  near  a  re-entrant 
comer. 

It  is  interesting  to  note  that  the  singular  nature  of  the  pressure  at  a  re-entrant  corner 
has  not  been  appreciated  by  computational  fluid  dynamists.  For  example,  Roache  [9] 
suggests  that  if  two  values  are  obtained  for  the  pressure  at  such  a  corner,  then  their 
difference  should  be  used  as  a  measure  of  truncation  error.  He  claims  that  a  non-single 
valued  pressure  is  non-physical.  Computationally,  for  the  Stoke  problem,  the  difference  in 
the  two  pressure  values  given  by  (3.3a)  and  (3.3b)  tends  to  be  unbounded  as  the  mesh  is 
refined  (see  Table  III). 

4.  Check  of  Accuracy 

As  we  said  in  §3,  the  finite  difference  method  we  used  to  solve  the  Navier-Stokes 
equations  is  second-order  accurate.  The  fact  that  the  pressure  possesses  a  singularity 
at  the  re-entrant  corner  may  affect  the  accuracy  of  the  scheme  near  the  step.  To  check 
the  accuracy  we  used  values  of  the  stress  function  at  the  wall,  and  at  different  positions 
downstream,  since  it  represents  a  sensitive  quantity  in  the  computation.  The  stress  at  the 


wall  is  defined  by 


dv 

7  =  ~P+  ' 

dy 


(4.1) 


where  v  is  the  transversal  velocity.  Let  ru  and  t 1  be  the  stresses  at  the  upper  and  lower 
wall  of  the  stepped  channel,  respectively,  (see  Figure  l)  then  in  order  to  eliminate  the 
pressure  constant,  we  define  the  value  At  by 


At  =  tu  -  t1  . 


(4.2) 


The  check  consists  of  taking  three  different  values  Ar^  ,  A r^,  ,  A Th3  corresponding  to 
the  computed  values  of  (4.2)  using  different  grid  sizes  hi,h2,h3  ,  where  h*  =  khk-i  and 
k  =  2,3  .  Then  we  assume  that  At  can  be  represented  as 

A rh  =  Ar  +  Chm  +  0(h*)  ,  (4.3) 


where  q  >  m  and  C  is  a  constant.  By  neglecting  smaller  order  terms  in  (4.3)  we  have  that 
the  order  m,  is  given  by 


m  =  log 


/  A Th2  -  A Thi  \  / 
V  A^3  -  AThJ  / 


log  2 


(4.4) 


We  made  the  computations  in  the  domain  with  the  inflow  boundary  at  x  —  —1  and 
the  outflow  boundary  at  x  =  1  ,  where  the  specification  of  Poiseuille  flow  was  used  at 
inlet  and  outlet.  The  test  positions  were  taken  at  i  =  ^,z  =  and  x  =  |.  We  used  two 
different  initial  grids  of  size  hi  =  1/12  and  hi  =  1/16,  respectively.  The  resulting  accuracy 


<r> 


is  given  in  Table  II. 
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TABLE  II 

Accuracy  for  the  Stokes’  problem  in  the  stepped  channel 


We  are  not  aware  of  any  previous  work  involving  the  computation  of  incompressible 
fluid  flow  in  domains  with  re-entrant  corners  in  which  a  similar  accuracy  check  had  been 
carried  out.  This  check  also  demonstrates  that  away  from  the  step  the  numerical  scheme 
retains  the  second  order  accuracy.  The  fact  that  the  accuracy  is  so  high  at  x  equal  to 
3/4  is  probably  due  to  the  simple  nature  of  v  and  p  in  Poiseuille  flow  and  also  due  to  the 
third-order  accurate  difference  formulas  for  the  divergence  of  the  velocity  and  the  pressure 
gradient. 


5.  Results 


x»-l  *=s 0  -  X sb  4 


Figure  1:  A  Channel  with  a  step. 

We  computed  flow  in  a  channel  and  in  an  axi-symmetric  pipe,  both  geometries  have 
re-entrant  corners  (Figure  1),  the  Reynolds  numbers  were  in  the  range  0  <  R  <  50 
we  use  the  pressure  corner  strategy  discussed  above.  The  results  include  graphics  of  the 
pressure  along  the  step,  (i.e.,  x  =  0,  0  <  y  <  1)  and  also  values  of  the  pressure  gap 
which  is  [p]  :=!  PCiX  -  Pc  y  |  .  All  the  computations  were  done  on  the  VAX  11/780  at  the 
Mathematics  Research  Center  at  the  University  of  Wisconsin-Madison.  To  check  the  self 
consistency  of  the  approach  we  introduce  the  following  device.  Since  near  the  corner  the 
pressure  behaves  like  p  ~  A^-456  +  »)  ,  with  A,  a  constant  depending  on  the 

angle  0  ,  then  the  gap  [pi  should  behave,  approximately,  as 

h°-456[p]=  C  j-  0(/i0'364)  ,  (5.!) 

with  C  a  constant  and  h  the  grid  size.  Table  III  shows  the  values  of  of  the  pressure  gap 

10 
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[p]  and  the  values  for  /i°'456[p]  ,  at  ten  different  grid  sizes  using  the  same  domain  and 
boundary  conditions  as  in  §4. 

TABLE  III 

Pressure  gap  for  different  grid  sizes 


h 

I?] 

h°-456ipj 

1 

12 

65.3894 

21.08119 

1 

16 

71.5178 

20.22504 

1 

20 

77.1634 

19.71243 

1 

24 

82.4458 

19.38327 

1 

28 

87.4137 

19.15760 

1 

32 

92.0902 

18.99143 

1 

36 

96.5233 

18.86576 

1 

40 

100.7433 

18.76781 

1 

44 

104.7752 

18.68960 

1 

48 

108.6414 

18.62613 

The  agreement  with  the  prediction  (5.1)  becomes  better  as  the  mesh  size  is  refined,  however 
since  the  lower  order  term  decays  very  slowly  the  third  column  in  Table  III  tends  slowly 
to  a  constant  for  very  fine  meshs. 

Figures  2. a  and  2.b  show  the  pressure  at  the  step,  R  =  0  ,  for  two  dimensional  flow 
in  a  channel  and  axi-symmetric  flow  in  a  pipe  using  grids  of  12  x  12,  24  x  24  and  48  x  48 
points,  respectively.  Note  that  for  the  pipe  flow  problem  the  Stokesian  pressure  (Figure 
2.b)  presents  the  same  type  of  behavior  as  for  the  channel  flow  problem,  i.e.,  that  the 


■»t  •«* 
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pressure  gap  becomes  unbounded  as  the  mesh  size  becomes  finer.  Figures  3.a  and  3.b 
show  the  pressure  at  the  step  for  R  greater  than  zero,  the  computations  were  performed  at 
Reynolds  numbers  equal  to  30,  40  and  50,  for  the  channel  flow  problem  and  at  Reynolds 
numbers  equal  to  10,  20  and  30  for  the  pipe  flow  problem.  From  Figures  3. a  and  3.b,  it 
can  be  observed  that  as  the  Reynolds  number  becomes  larger  the  pressure  gap  tends  to 
be  smaller.  It  seems  to  be  that  the  assumption  that  near  the  corner  the  solution  of  the 
Navier-Stokes  equations  can  be  expanded  in  terms  of  Moffatt’s  eigenfunctions,  is  only  valid 
for  small  Reynolds  numbers. 

6.  Conclusion 

We  have  presented  a  method  of  treating  the  pressure  at  a  re-entrant  corner,  this 
strategy  is  motivated  by  the  need  for  obtaining  pressure  values  at  the  wall.  The  resulting 
pressure  value  at  the  corner  may  be  regarded  as  double  valued.  These  values  of  the  pressure 
at  the  re-entrant  corner  are  used  to  evaluate  the  pressure  gradient  at  adjacent  grid  points. 
This  strategy  gives  results  which  are  consistent  with  decreasing  mesh  length  and  also  can  be 
used  to  treat  axi-symmetric  and  three-dimensional  viscous  corner  flow  without  introducing 
any  modification  in  the  nature  of  the  algorithm.  We  used  as  test  problems  channels  and 
axi-symmetric  geometries  with  re-entrant  corners  and  computations  were  performed  at  low 
and  moderate  Reynolds  numbers.  This  way  to  treat  the  pressure  at  a  re-entrant  corner 
seems  to  work  the  best,  among  a  few  other  strategies  we  used,  applied  to  the  numerical 
scheme  in  [10].  The  numerical  results  show  that  the  pressure  singularity  can  be  modeled 
very  accurately,  making  the  computation  more  efficient. 
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